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Standard algorithms for the numerical integration of the Langevin equation require that inter- 
actions are slowly varying during to the integration time-step. This in not the case for hard-body 
systems, where there is no clear-cut between the correlation time of the noise and the time-scale 
of the interactions. Starting from a short time approximation of the Smoluchowsky equation, we 
introduce an algorithm for the simulation of the over-damped Brownian dynamics of polydisperse 
hard-spheres in absence of hydrodynamics interactions and briefly discuss the extension to the case 
of external drifts. 
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I. INTRODUCTION 

The discovery that suspensions of colloidal particles can be tuned to be excellent experimental realizations of the 
idealised hard-sphere (HS) system (l6l . [26[ has triggered in the last decade a renewed interest in the theory and 
simulation of hard-spheres. Since colloidal hard spheres have a radius in the size range of 10 ~ 1000 nm, they can 
be at the same time much bigger than the solute particles and small enough to have enough thermal energy to 
disregard gravitation; hence colloidal HS dynamics can be modelled as Brownian motion in presence of hydrodynamic 
interactions. From the theoretical point of view, hydrodynamics force are often disregarded and the simple model 
of Brownian HSs is employed to understand real HS suspensions; yet, even for the simple Brownian model only 
approximate theories are possible and simulations are needed to discriminate among them. 

Brownian dynamics (BD) algorithms integrate numerically Langevin equations; a common requirement of such 
algorithms is that interactions in the system should vary little during an integration time-step. Under such assumption, 
particle displacements are calculated keeping forces constant during the integration time-step 0. In the case of 
continuous potentials, computational efficiency worsens as the interactions become steeper. In the extreme case of 
hard-body interactions, stochastic calculus is not naively applicable and standard numerical integrators become ill 
defined. On the other hand, the Kramer's equation [Tl| associated with Brownian motion is well defined when suitable 
boundary conditions taking account of stepwise interaction are implemented (l9| . One general strategy to develop 
numerical integrators for stochastic differential equations is to work on the short time expansion (in particular on 
Trotter expansions [25]) of the associated Fokker-Plank equation Q. We will develop an approach similar in spirit 
[20^ . by coming to integrate the over-damped Brownian dynamics (OBD) with stepwise interactions via suitable 
approximations of the associated Smoluchowsky equation (SE). 

In section |H] we recapitulate and justify the standard event driven BD algorithms for homogeneous systems of HSs; 
in section Hill we extend such schemes to the case of polydispersity and in section HVl we analyse the extension to the 
case of constant drifts. 



II. OVERDAMPED BROWNIAN DYNAMICS 

The stochastic differential equation describing an homogeneous system of overdamped Brownian particles is 

d t n = fi + ii 

where 7^ are the coordinates of the i th particle, fi are the non-random forces acting on i and £j is an uncorrelated 
Gaussian noise 

(ii ®£j) = 2D5 af} 

V / a/3 

whose amplitude is twice its diffusion coefficient D. In the case of HSs of diameter a, the force fi contains infinite 
impulsive contributions due to the HS interaction potential 

y = f for \ fy \ > a 
13 1 00 for \ fij I < (7 

Such a force is not Lipschitz continuous and standard methods for stochastic differential equations become out of 
reach [H, [HJ. Notice that already for systems of classical particles of mass m in the micro-canonical ensemble, the HS 
force take a peculiar velocity-dependent form 

fijdt = —mVijS (\f i: j \ - a) 

where v\j is the relative velocity between particles i and j along the direction of = r j — fj ; in the case of overdamped 
Brownian HSs, velocities are not defined and the HS conditions \ fij \ > a must be interpreted as boundary conditions. 
In fact, the Fokker-Plank equation associated to the OBD of HSs takes the very simple form of a free Smoluchowsky 
equation [22] 

d t P(v,t)=DV 2 P{r,t) (1) 

with suitable boundary conditions; here P (r, t) is the probability distribution function (PDF) for the positions 
r = {fi}. It is an equation of the form of a divergence <9 t P = div (j) in the current j = Dd r P with div (j) = d r ■ j; all 
the complexity is in the implementation of the hard-sphere impenetrability by a reflecting (zero current) condition 
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on the time-dependant boundary d£l corresponding to |fy (t)\ — a (i.e. spheres i and j are at contact at time t); here 
and n is the normal to d£l. 

To build up an algorithm to integrate such a system, one has to rely on physical intuition: considering integration 
steps At small enough, particles will perform on average free random walks until some couples of particles are "near 
enough" to interact. This is the basis of many algorithms for OBD in the case of HSs: first, independent particles 
displacements are extracted according to the free Green's function for single particle diffusion 

G{ ree (r,t + At\f ,t) oc exp [- (r - f f /2DAt\ (2) 

; then, overlaps are taken account to correct such displacements 0, 0, @, H3, HH HH H3|- In ah such schemes, the 
implicit assumption is that for small time-steps At the evolution of the full P (r, t) factorizes either in single particle 
free evolutions p (rj, i) or in the evolution p (fi, fi 7 t) of two interacting particles. 

As shown in [231. naively chosen corrections can lead to the wrong dynamics. For purely HS interactions, a naive 
algorithm [|| |24| that transforms the displacements Ari in Active velocities V{ = Afi / At and evolves the system 
according to the rule of standard event-driven molecular dynamics [17] (EDMD) has been shown to approximate 
correctly the SE of the system [2(j. In such an approach, the time step At is fixed; at each time-step, the velocities 
of the particles are extracted according to the Maxwell distribution at a Active temperature T and a fully fledged 
EDMD simulation [l7| is performed between time t and t + At. The temperature T is chosen such that the average 
displacement in absence of collisions is exactly eqf2J 

In order to justify such algorithms, several hypothesis must be done. First, the time step At must be small enough 
that only binary collisions must be relevant, i.e. the average displacement { | Ar^ | } ~ At 1 / 2 must be much smaller than 
the average inter-particle distance 

(\An\) ^p- l ' d -a (3) 

; here p is the number density and d is the dimension of the system. In such a limit, the interaction among two 
overdamped Brownian HSs i and j can be mapped to the problem of a point overdampcd Brownian particle in 
presence of a sphere by the change of coordinates (fig. [TJ 

Tij = n- fj 

Rcm - + { ' 
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; in such a reference system, the Brownian center of mass Rcm is subject to free diffusion d t RcM — S while 
satisfies the SE with spherical boundary conditions 

dtfn = e (5) 



|f| > a 

; here £ = £j — £j and 3 = k* + /2- Notice that £y and S are mutually orthogonal Gaussian noises and therefore 

the equations for Rcm and for fy can be solved independently. Equation ([5]) can be exactly solved P, 0| but the 
solution is in the form of an infinite sum in the Laplace domain; therefore, further approximations are needed as it 
is not suitable for the fast for numerical implementations necessary to simulate many-body system. In particular, 
the condition of small displacements during the time step At can be pushed to satisfy also an additional "flat wall" 
condition 

(|Ar-|)«a (6) 

In such a situation, binary collisions modelled by eq.([5]) happen on average only between particles at an initial 
distance \ fij\ = a, i.e. at distances much smaller than the radius of curvature of the spherical boundary. In such a 
situation the boundary can be approximated as a flat wall. Shifting to Cartesian coordinates such that the origin lies 
on the intersection of fij with boundary and orienting the y.z axis tangentially to the surface, the system factorizes 
in two free Smoluchowsky equations for the y,z coordinates and a one dimensional equations for the 
coordinate 




which can be exactly solved [22] with the image method [18] (figEJ) 



-l-r-i> l 2 -|P+Pq| 2 

11 (f) re < e 2DAt + e 2da* for x < V 



for x > 

The whole solution G™ aU for a point particle starting in r*o consists of the superposition in the x < semi-space 
of the free Green's function ([2]) of a particle in fo = (xq, yo, Zq) and an image particle in Tq = (— xq, yo, zq). Such 
a solution can be implemented with just a single operation by extracting the new position r(t + At) according to 
(J2) and reflecting the x coordinate whenever x > 0. An event-driven algorithm implementing such scheme is the 
following: 

1. extract the random displacements Ax, Ay, Az, 

2. define a Active velocity v x = Ax/ At 

3. calculate the fictive collision time t c : xq + v x t c = 

4. calculate fictive the post-collision velocity v* = —v x 

5. calculate the final position x (t + At) = x (t) + v x ■ t c + v* ■ (At — t c ) 

When mapping back from the Brownian Center of Mass (BCoM) reference system to the original particles' coordi- 
nates, the whole Brownian collision fi (t) — > fi (i + At), rj (t) — > fj (t + At) follows a procedure strictly recollecting 
the collision of two classical HSs: 

1. extract two random displacements Ar\-, Afj according to ^ 

2. define two fictive velocities Vi = Af^/At, Vj = Arj/At 

3. calculate the fictive collision time t c 6 [0, At], and the normal n* between the two spheres at contact at time t c 

4. calculate the fictive post-collision velocities v* = ui — 2 (n* ■ Vij ) n* and v* — Vj + 2 (n* • t% ) h* with Hij = vi — Vj 

5. calculate the final positions fi (t + At) = fi (t)+Vi -t c +v* ■ (At — t c ) and fj (t + At) = fj (t)+Vj -t c + v*- (At — t c ) 

Therefore, any event-driven molecular dynamic code [171 ] can be adapted to simulate the OBD of HSs by extracting at 
each step At the velocities of the particles according to independent Gaussian distributions such that (yf ) = 2dDAt~ 1 \ 
a consistency check to perform is to ensure that within the chosen time-step At less than one collision per particle 
occurs on average. 
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Reflecting Wall 

Figure 2. Image method of solution for the Smoluchowsky equation in presence of a fiat reflecting boundary. The probability 
of an overdamped Brownian particle of reaching a point f at time t is given by the sum of the probability G[ ree (r, t + Ai|r*o, t) 
that the particle goes from fo to r by free diffusion plus the probability G{ ree (r,t + At\rg,t) that the particle goes from the 
image of the initial point to the same f also by free diffusion. The latter is equal to the probability Q^ ree (ff* ) t + At\fo, t) 
that the particle goes by free diffusion from rb to the image r* of the point r. Therefore, to implement numerically the image 
method, it suffices to implement the following algorithm: (I) extract the final position r according to the solution G[ re£ of the 
free Smoluchowsky equation and (2) reflect f if it goes beyond the hard boundary. 



III. POLYDISPERSITY AND BROWNIAN COLLISIONS 

In realistic colloidal system there is an inherent polydispersity in particle size; more generally, one could be also 
interested to mixtures of HSs with different characteristics as in all the studies where crystallization must be avoided. 
Let us now suppose that our system is composed of HSs of diameters <jj subject to overdamped Brownian motion 
with particle-dependent free diffusion coefficient Di, i.e. to noises of amplitude (|^) = 2dDi, as usual, trajectories 
are subject to no-flux boundary constraints \fi(t) — ^j(t)\ > Cy where cry = (<7j + <jj) /2. Let's suppose again to fix 
a time-step At small enough to consider only binary collisions; for two particles i and j the equations become 

d t n = & 

\nj \ > <Jij 

The first step is to separate equations ([7]l: transformation (U) does not succeed any longer as the transformed 
noises have a non-zero correlation Di — Dj. To properly define the "Brownian center of mass" (BCoM), we start from 
the ansatz Rcm — a i^i + a j? and impose zero correlation among the random displacements of the BCoM and the 
inter-particle distance: = ^AiZcjif Ar^\ oc aiDi — djDf, a proper dimensionless choice is <Zj = (D{ + Dj)/ Dj cc D^ 1 . 

In the limit of At small enough such that the jr^ | = cry the boundary can be approximated with a fiat hard wall and 
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the zero flux problem can be solved along the dimension perpendicular to the wall (we assume it is the x direction). 
Let's again define Active velocities both in the original reference system and in the BCoM system v a — Ax a / At for 
a e {i,j,CM,ij}. 

The change of coordinates to the BCoM system is 



VCM 



= A 



Di+Dj Di+Dj 



I 



-1 



with inverse 



= A~ 



VCM 



' D,+D, 
Di 



VCM 



(Di+Dj) 2 D % +D 3 

In the BOcM system the collision corresponds to imposing the no-flux boundary collision and is simply 



V CM 



1 

-1 



and therefore the fictive velocities after the collision are 



with collision matrix 



A- 1 



J CM 



= A- 1 



I 

—I 



VCM 



VCM 



A- 1 



1 

-1 



A 



Vi 



C 



Brown 



1 

—I 



A 



Di-Di 



2Dj 



Di+Dj Di + Dj 
2D, Di-Di 



Di+Dj Di + Dj 

This is to be compared with the classical collision matrix for two elastic particles of masses m i; mj 



Cclass — 



2m, 



mi+rrtj rrti+mj 
2m; m,-roj 



rrii+rrij rrii+rrij 



that has a similar structure if one but with switched indexes i <H- j such the role of the mass m; during a collision 
is played by the inverse diffusivity £)r . 

As a check, we consider a fixed particle with Di = 0; this is equivalent for particle i as having an infinite mass 
(noise does not move it) and the collision is 



C 



1 

2 —I 



Vi 

2vi — Vi 



as it should be (notice that v* = —vj as Vi = for Di = 0); an analogous result comes by sending the "mass" of 
particle j to zero (i.e. Dj = oo). 

To summarize, let's recall that the full event-driven collision scheme for classical particles is 



1. calculate collision time t c from the "good" root of \\fij + Vijt c \ 

2. bring particles at contact f. \ = r, L + Vit c ,fj = fj + Vjt c 



3. let Sij = fij(t c ) , <Jij = \\(Tij\\ ,Vi = Vi ■ &ij,Vj 



4. pre-collision: 



Vi 



a 



class 



5. collision: v[ = Vi - Vi&ij + v-frij , Vj = Vj - Vj&ij + v'j&ij 



Therefore, to modify an Event Driven code for polydisperse HSs into an Event Driven Brownian Dynamics, simply 
extract random displacements Afi at fixed intervals t, t + At,t + 2At, . . ., define fictive velocities Vi = Afi/ At and 
evolve the system for a time At using for the collision the matrix C ' Brown instead of C c i ass . 
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IV. CONSTANT DRIFTS 



Insofar, only system not subject to external forces have been considered. For overdamped Brownian motion, 
constant forces add constant drifts to the random displacements of the particles. Polydisperse particles are to be 
expected to experience drifts of different magnitudes even in presence of an homogeneous fields (like gravity or an 
electrical field). Therefore, the two body equation becomes 



f dt?ij = 4 + 9ij (8) 

where = gi — gj is the difference among the constant drifts of the two particles. Notice that homogeneous drifts do 
not produce any change in the equations for the inter-particles distances as g^ = 0, but just add a constant drift to 
the BCoM Rcm] therefore, Brownian Event Driven simulations can be implemented extracting particle displacements 
according to the Green function of the Smoluchowsky equation for a single particle with drift instead of @. 

To see the effects of no-zero g^-, let's consider again a At small enough such that the collision can be approximatively 
by a flat wall. Factorizing the motion in the directions perpendicular (x) and parallel (y, z) to the wall, one is left 
with solving the Smoluchowsky equation with a reflecting boundary in the case of constant drift. This problem has 
been solved at the beginning of the last century in the seminal paper by M. Smoluchowsky [22| (pp. 569-574; see [l4j . 
pp. 2714 for an English version). Assuming that the reflecting boundary is the plane dil = {x = 0}, the evolution of 
the probability distribution function follows the equation 



f d t p = Ddlp + cd x p 

\ subject to (d x +c)p(x,t)\ x=0 

where c = — /3g, j3 — 1/fceT is the inverse temperature and g is a constant force. 
Solving for the initial condition p(x,t = 0\xq) = 8 (x — xo) one obtains the solution 



p(x,t\x ) 



2VrrDt 



4Dt -f e 



(*+*or 



2D 4D 



(9) 



where er/c (z) — I— erf (z) is the complementary error function and the error function is erf(z) = \J A/n J* ds exp [— s 2 ] . 

Such a solution is not amenable of a simple geometric implementation in terms of a naive collision mechanism: in this 
case, for each collision one must transform the coordinates to the BCoM reference system, extract the displacements 
according to ^ and transform back to the original coordinate system. As a further caveat, high enough constant inter- 
particle drifts imply the accumulation of particles at short distances (fig El); the appearance of such inhomogeneous 
structures is a critical situation is critical for event-driven algorithms as it can produce unacceptable slowing-down of 
the simulation (via the growth of the number of collision per unit time) and eventually numerical errors [17fl • Notice 
that such issues of non-zero drifts among nearby HSs are often disregarded in the simulations of sheared particles 
where collisions between hard disks or hard spheres are implemented as elastic collisions of the fictive velocities 

a Eta Em . 

Since for HSs structural quantities like the pressure are strictly related to the radial distribution 
function at contact, a careful analysis of the importance of the drift term in relation to the strength of the noise 
should be performed to avoid disregarding possible relevant corrections. 
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